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ABSTRACT 

We investigate the effects of viscosity on disk-planet interaction and discuss 
how type I migration of planets is modified. We have performed a linear calcula- 
tion using shearing-sheet approximation and obtained the detailed, high resolu- 
tion density structure around the planet embedded in a viscous disk with a wide 
range of viscous coefficients. We use a time-dependent formahsm that is useful in 
investigating the effects of various physical processes on disk-planet interaction. 
We find that the density structure in the vicinity of the planet is modified and 
the main contribution to the torque comes from this region, in contrast to invis- 
cid case. Although it is not possible to derive total torque acting on the planet 
within the shearing-sheet approximation, the one-sided torque can be very dif- 
ferent from the inviscid case, depending on the Reynolds number. This effect 
has been neglected so far but our results indicate that the interaction between 
a viscous disk and a planet can be qualitatively different from an inviscid case 
and the details of the density structure in the vicinity of the planet is critically 
important. 

Subject headings: planet and satellites: formation — solar system: formation 
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1. Introduction 

Disk-planet interaction and type I planetary migration is one of the main issues in the 
theory of planet formation. A low mass planet or the core of the gas giant planet embedded 
in the protoplanetary disk excites spiral density wave in the disk, which carries angular 
momentum flux away from the planet. The linear calculation of the spiral density wave is 
first carried out by Goldreich and Tremaine (1979), and there has been a number of extensive 
work assuming inviscid, isothermal disk such as Ward (1986, 1989), Artymowicz (1993), and 
Tanaka et al. (2002). In the case of an inviscid, isothermal disk, protoplanets migrate 
towards the central star as a result of disk-planet interaction very rapidly, thereby imposing 
a serious problem on the theory of planet formation. 

The decay of the orbital semi-major axis of the protoplanet is a result of the subtle 

difference between the inner disk (the part of the disk inside the orbit of the protoplanet) 
and the outer disk (the part of the disk outside the orbit of the protoplanet). The outer 
disk exerts torque on the protoplanet in such a way that the semi-major axis decreases, 
and the inner disk exerts opposite torque. In the lowest order of the local approximation, or 
shearing-sheet approximation, these two effects cancel each other, and the difference appears 
at higher order. Therefore, type I migration can be very sensitive to the physical state of the 
disk. Indeed, when non-isothermal effect is taken into account, the direction of migration 
can be reversed (Paardekooper and Mellema 2006, Baruteau and Masset 2008, Paardekooper 
and Papaloizou 2008) . It is also discussed that the global magnetic field can slow down the 
rate of type I migration or even reverse the direction (Terquem 2003, Promang et al. 2005, 
Muto et al. 2008). 

In order to construct a predictable, self-consistent model for planet formation, it is 
important to make clear how different physical processes of the disk affect disk-planet inter- 
action and type I migration. In this paper, we revisit the effects of viscosity on disk-planet 
interaction. There has been a number of work on this topic. Meyer-Vernet and Sicardy 
(1987) discussed that the viscosity does not affect the torque exerted on the planet as long 
as its value is small, although the shape of the wake may be modified. Papaloizou and Lin 
(1984, and a series of their paper) and Takeuchi et al. (1996) calculated that the damp- 
ing of the density wave and the transfer of the angular momentum from the density wave 
to the background disk, resulting in the gap formation around a massive planet. Masset 
(2001) investigated the effects of viscosity on the horseshoe regions and discussed how the 
corotation torque may be altered by the effect of viscosity. More recently, Paardekooper and 
Papaloizou (2009b) investigated the dynamics of corotation region using both linear and 
non-linear calculations. 

Although there has been a number of investigations on the disk-planet interaction in 
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a viscous disk (for example, Masset 2001, 2002 or Paardekooper and Papaloizou 2009a for 
corotation torque, D'Angelo et al. 2002, 2003 for high resolution numerical study), there 
has not yet been a study of wide range in the viscous coefficient that requires an analysis 
on the detailed density structure in the vicinity of the planet. In this paper, as a first step 
for the complete investigation, we show results of linear calculation in a local, shearing-sheet 
analysis. We use time-dependent methods to calculate the disk response against the planet 
potential. Although this method has not been widely used so far, this is useful in investigating 
the effects of various physical processes on disk-planet interaction. We calculate the non- 
axisymmetric density structure around the planet and investigate how the resulting torque 
is altered by the effect of shear viscosity. We have studied wide range of the parameters 
of viscous coefficient and calculate the density structure with high resolution. We find that 
the density structure in the vicinity of the planet is altered in a viscous disk, with viscous 
coefficient of ~ 0.1 in terms of a. (standard a parameter, see equation (1271) for definition), 
which may be realized as a result of the turbulence induced by magneto-rotational instability 
(MRI, Balbus and Hawley 1991, Sano et al. 2004). 

Since it is not possible to calculate the total torque (differential torque between the 
disk inside the orbit of the planet and outside) in shearing-sheet approximation, we look at 
one-sided torque to investigate how type I migration rate can be affected qualitatively by 
viscosity. One-sided torque is the torque exerted by one side of the disk with respect to the 
planet's orbital radius. We also note that shearing-sheet calculation is a useful first step 
before differential torque is calculated (Tanaka et al. 2002). We find that as we increase 
the values of viscous coefficient a, one-sided torque stays unchanged until a <^ 0.01, then 
increases until a ~ 1, and finally decreases. The torque can be factor of two larger than 
inviscid case when a = 0.1, and more importantly, the enhancement of the torque is a 
result of the modified density structure in the vicinity of the planet, which has not yet been 
investigated in detail. Our results indicate that the physical mechanisms of the disk-planet 
interaction in a viscous disk may depend on the detailed density structure around the planet. 

The plan of this paper is as follows. In section [2] we show the basic equations and 
describe the time-dependent methods to obtain a stationary, non-axisymmetric pattern of 
linear perturbation. We also discuss how two-dimensional (2D) and three-dimensional (3D) 
analyses are related in this section. In section[3l we show the results of 2D and 3D calculations 
and show how the torque behaves as we vary the amount of viscosity. In section HI we 
explain the qualitative behavior of the torque using a simple analytic model. Section [5] is for 
discussion and summary. 
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2. Basic Equations and Numerical Methods 

2.1. Linear Perturbation Analysis 

We consider isothermal Navier-Stokes equations with one planet using shearing-sheet 
approximation. 

^^+V-{pv) = (1) 
dv 1 

— + v-\/v + 2npe, xv = Vp + Snlxe^ + uV^v + -z/V(V ■ v) - Vtpp (2) 

where p is density, v is velocity, Qp is the Kepler angular velocity of the planet, ipp is the 
gravitational potential of the planet, u is shear viscosity. The third term of left hand side 
of equation ([2]) is Coriolis force and the second term of right hand side of equation ([2]) is 
tidal force. We assume the planet is located at the origin and stationary with respect to this 
coordinate system, i.e., 

il^p = tppix, y,z), (3) 

where ipp does not depend on time. In this paper, we neglect bulk viscosity for simplicity. 
We also neglect vertical stratification and assume that the background density without a 
planet is homogeneous. This gives an uncertainty in the box size in the 2;-direction, which 
will be discussed in Section 12. 2[ We also neglect the effect of global gas pressure gradient 
exerted on the background disk and assume that the background gas is rotating at Kepler 
velocity. 

The unperturbed state without a planet is given by 

p = Po = const (4) 
3 

vo = --ilpxey. (5) 

In the presence of viscosity, if we perform a global analysis, there is a mass accretion onto 
the central star in general. However, in the shearing-sheet approximation, where linear 
background shear is assumed, this effect is not taken into account. We expect that the 
density structure only in the vicinity of the planet can be well approximated even if we 
neglect the global mass accretion. 

We consider linear perturbation. All the perturbation quantities are denoted with S, 
e.g., 6p for density perturbation. Perturbation equations are given by 
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dt 2 ^ dy J p y ^ 2 P ^ ^ 

= -c^V^ + uV^5v + -z/V (V • 5v) - Vtpp (7) 
Po 3 

We solve equations (Q and ([7j) to obtain a steady state solution and calculate torque exerted 
on one side (either a; > or a; < 0) of the disk by the planet. The torque exerted on the 
planet is obtained as a backreaction of this torque. Since we use shearing-sheet approxi- 
mation, torque exerted from each side of the planet is the same in magnitude and opposite 
in sign. Although we do not obtain a net torque, it is still possible to investigate how the 
disk structure is affected by the viscosity and qualitatively predict how the disk and planet 
interact. 

In order to obtain a steady state solution, we solve linear perturbation equations ([6]) and 
([7]) using Fourier transform methods given by Goodman and Rafikov (2001). We transform 
the equations into the shearing coordinate (t', x', y', z') defined by 

t' = t (8) 
x' = x (9) 

y' = y + hlpxt (10) 

z' = z. (11) 

In this coordinate system, temporal and spatial derivatives are given by 

dt dt' 2 " dy'' ^ ' 

di~d^'^2''^d^" ^^^^ 



d d d d 
dy dy'' dz dz' 
Then, perturbation equations become 



(14) 



|^^ + V5t; = (15) 
of Po 



d 1 

— 5v - 2VLp5vyen, + -VLpSv^By 

-c^V— + pW^5v + l-uW (V ■ 5v) - V^p, (16) 
Po 3 
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where the spatial derivatives in V are given by equations f|T3|) and f|T^ . 

The coefficients of equation ([T^ and ffTBl) are now independent of (x', y' , z'). Therefore, 
if we Fourier transform in spatial directions, we obtain a set of ordinary differential equations 
decoupled for each (A;^, fc^, k'^) mode (Goldreich and Lynden-Bell 1965), 

5f{t\ x', y\ z') = J2 Sfit'^ K) exp [tiKx' + k'^y' + k'^z')] . (17) 

The relationship of wavenumber in (x, y, z)-coordinate and (x', y' , z')-coordinate is given by 

k^it) = k',+ hl^kyt. (18) 

~ — (19) 

Equation fllSp indicates that the value of radial wavenumber in (x, z)-coordinate evolves 
with time owing to the background shear. The value of radial wavenumber in shearing 
coordinate k'^ gives the initial value of radial wavenumber in (x, z) -plane. Equations of 
continuity ([6]) and motion ([7]) are now 

d S p 

— h ikx{t)5vx + ikydvy + ik.^5vz = 0, (20) 



^Svx - 2^lp6vy = -c^ik^{t)— - v{k^{tf + kl + k'l)5v^ 
dt pq 

1 

-i:J^kx{t){kx{t)6vx + ky6vy + k^^v^) - ikx{t)ipp, 



(21) 



^^^y + ^^P^'^x = -c^ikyy- - z/(A;^(t)^ + + ^l)^'^y 
-\i>ky{k^{t)5v^ + kydvy + k:,5vz) - ikytpp, (22) 



and 

^Sv, = -cHk,^ - v{k^{tf + kl + kl)5vy 
dt po 



--uky{k^{t)6v^ + ky5vy + kz5vz) - ik^tpp. (23) 

In equations ( !20l) - (j23l) . we write all the terms in (t, x, 2;)-coordinate. These equations 
describe the excitation and evolution of density wave by the source ipp- The source of the 
density wave becomes zero when wavenumber {k'^, k'y,k'^) is large. 
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As an initial condition, we assume that there is no perturbation: 5p{t = {]) = 5v{t = 
0) = 0, and k'^ is taken to be sufficiently large in absolute magnitude and negative (positive) 
in sign for positive (negative) ky. If we use spatial resolution of x-direction Ax, we should 
take k'^ = =F27r/Ax, where upper (lower) sign is used for positive (negative) ky. As a result of 
time evolution, k^it) increases (decreases) for each positive (negative) ky mode. When kx(t) 
reaches ±27r/Ax, where upper (lower) sign is for positive (negative) ky modes, we obtain a 
steady state profile of perturbation quantities in Fourier space for non-axisymmetric {ky ^ 0) 
modes. The Fourier amplitude of each k^ mode in steady state is given by the time evolution 
data through equation (|T8l) . The profile of physical quantities in Fourier space is then inverse 
Fourier transformed to obtain values in real space. Since we are interested in torque, we do 
not calculate axisymmetric modes iky = 0). 

The standard procedure to obtain steady state solution is as follows. Firstly, stationary 
solution in the frame corotating with the planet is assumed: d/dt = 0. Then, Fourier 
transformation in the y- and 2;-directions is performed to obtain the ordinary differential 
equations in the x-direction. Finally, these ordinary differential equations are solved by 
imposing outgoing boundary condition, or equivalently the boundary condition that admits 
only trailing wave, see e.g., Goldreich and Tremaine (1979), Korycansky and Pollack (1993), 
or Tanaka et al. (2002). 

In the presence of viscosity, this method introduces higher order derivative with respect 
to X, since viscous terms include the second order derivative in the radial direction. The 
resulting ordinary differential equations are higher order in the x-derivative than in inviscid 
cases. The terms with the highest order derivative comes from viscous terms and their 
coefficients are viscous coefficient u. In this case, it is difficult to take natural — > limit, 
since equations become singular in this limit. 

Time- dependent approach overcomes this difficulty since the order of time derivative is 
not affected by viscous terms since they do not include any time derivative. It is also easy 
to take z/ — > limit since viscous terms simply drop from equations (l2Tl) - fl23l) in this limit 
and they do not introduce any singularity. 

In our formulation using the Fourier transform in sheared coordinate, the outgoing 



^ The reason why we obtain a steady state as a result of time evolution is as follows. Let us assume 
that if \kx\ > Kc, the source '0p is so small that we can approximate it to be zero. The Fourier amplitude 
of the specific mode in it, x, y, z)-coordinate at time t = to is given by the solution at t = to with 
k'x = — {3/2)npkyto mode. Therefore, the absolute magnitude of k'^ that contributes to the Fourier 
amplitude of fixed k^o becomes large as time to increases. If the absolute magnitude of k'^ is larger than Kc, 
there is no wave excitation for this mode until |fcx(i)| becomes smaller than Kf.. Therefore, when |/c^| > K^, 
the contribution to the k^Q mode is always the same regardless of t^. 
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boundary condition in the radial direction that is assumed in the standard stationary for- 
mulation is not exactly satisfied, since Fourier transform introduces the periodic boundary 
condition in the sheared coordinate. However, we find that outgoing boundary condition 
is effectively satisfied if we vary the box size in the x-direction depending on the modes 
specified by ky and kz. Details of our methods are given below. 

We also note that if we assume an initial condition with non-zero perturbation, addi- 
tional homogeneous wave is introduced. This wave has both leading {k^ky < 0) and trailing 
{k^ky > 0) components and the resulting solution is simply the superposition of the specific 
solution of Equations fl20|) - fl23|) assuming zero initial condition and homogeneous [ipp = 0) 
solution assuming the specified initial condition. The solution with non-zero initial condition 
is irrelevant in the present problem because we consider the perturbation that is induced by 
the gravitational perturbation by the planet. 

We write the box size in (x, ?/, 2;)-directions by {L^, Ly, Lz) and the coordinate system 
extends —Lx/2 < x < L^jl and so forth for other directions. Practically, our procedure to 
obtain the stationary, non-axisymmetric structure of the disk is summarized as follows. 

1. We take k'^ = ^2n/Ax for positive (negative) ky modes. 

2. Equations (!20|) - (l23|) are solved with initial condition 6p{t = 0) = 6v{t = 0) = 0. 

3. Resulting solutions in {k^, ky, k^) space are inverse Fourier transformed to real space. 

In steps 2 and 3, there are some points we need to take care. In calculating the Fourier 
modes, we have varied the mesh number in kx directions, which corresponds to the step 
size of time, in such a way that all the oscillations are well resolved for each {ky, kz) mode. 
The mesh number in k^ direction becomes as much as 10^ for small ky. We perform inverse 
Fourier transform in k^ direction using the increased number of mesh in k^ direction. This 
procedure effectively makes the box size of x-direction longer for a specific ky mode, and 
the data corresponding to necessary x, which is —Lx/2 < x < Lx/2, is then interpolated 
from the resulting profile in {x, ky, kz) space. This data is then used to perform inverse 
Fourier transform in {ky, kz) directions. In this way, it is possible to avoid aliasing effect 
in x-direction and outgoing boundary condition in the steady state solution is effectively 
satisfied. We do not use any window function, which has been incorporated by Goodman 
and Rafikov (2001), in Fourier transform. We have checked that the resulting torque is not 
affected by the window function in the absence of viscosity. Also, since window function 
introduces additional artificial dissipation, this affects our results with viscosity. 

Once we have obtained the profiles of perturbation quantities in real space, we can 
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calculate the one-sided torque exerted on x > part of the disk by the planet by 

fLx fLy/2 /•Lz/2 
-Ly/2 J-Lz/2 



r-Lx rLy/2 />Lz/2 gi 

T^-Tp / / dxdydz5p{x,y,z)^, (24) 

Jo J -Ly/2 J-Lz/2 



where is the semi-major axis of the planet, and the torque exerted on the planet is 
obtained as backrcaction. For later convenience, we define "torque distribution" T[x,y,z) 
that is simply the torque exerted on the fluid element located at {x, y, z), 

T{x, y, z) = -rp6p{x, y, z)-^ (25) 

and "torque density" that is the torque exerted on an annulus of the disk at x, 



T{x) ^-Tp I I dzdy5p{x, y, z)^ (26) 



dy 



We normalize equations using c, r2p, and po so that the homogeneous equations ('^p=0) 
contain only one dimensionless variable a that is defined by 

« = ^. (27) 

where H = c/Qp is the scale height of the disk. Since we consider the linear perturbation 
analysis, the amplitude of the perturbation is proportional to the normalized planet mass: 

p = GMp/Hc^. (28) 

We show the results with /x = 1 in subsequent sections. 

We use the fifth order Runge-Kutta method to solve time evolution and FFT routine 
given by Press et al. (1992) to perform Fourier transform. We use box size of = lOH 
and Ly = AOH, and the grid number in x- and |/-directions {Nx,Ny) = (512,512). The 
three-dimensional results depend on the box size in 2;-direction, L^, which is discussed in the 
next subsection. 



2.2. Vertical Box Size and the Difference between 2D and 3D Calculation 

In this subsection, we discuss how 2D and 3D calculations differ from each other and de- 
termine the relevant box size in the ^-direction with which the effect of vertical stratification 
is effectively taken into account. 
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We compare 2D calculation and 2D mode of 3D calculation. By "2D calculation", we 
mean the gravitational potential of the planet is given by 

^p,2d(x,i/) = -^=£%=, (29) 

where €20 denotes the softening length of 2D calculation, which is usually incorporated in 
most of 2D work. By "2D mode of 3D calculation" , we mean that the analysis is restricted 
to fc^ = mode of 3D calculation. Three-dimensional potential is given by 

V-p.sDlx, y, z) = =======. (30) 

+ y^ + z^ + 

Therefore, gravitational potential for "2D mode of 3D calculation" is given by 

1 /"^^/^ , GM„ 



'ipp,3B{x,y) = -— / dz 

J-'z J-LJ2 



'P 



L,/2 + ?/2 + 2;2 + 



GM 

— T^log 



(31) 

+ y^ + e^j3 

where we denote softening length in 3D calculation by csd- We note that ipp,2T) and '?/'p,3D 
coincides if + + ^ L^, but they differ if + + ^ L^, since in this case, 
3D potential behaves as V'p^sD ~ logr while 2D potential behaves as '?/'p,2D ~ '^l'^-, where 
= + + e^. Vertical averaging makes the potential weaker in the vicinity of the planet, 
thereby producing small perturbation in 2D mode of 3D calculation. Figure [1] shows the 
torque obtained by 2D inviscid {v = 0) calculation using 'ipp^2D and ipp^sB with various L^. 
We note that in producing Figure [H one-sided torque is calculated according to equation 
(121]) and the value is normalized by 

= fi'^poL^rpHc^. (32) 

It is clear that 2D and 3D results coincide if is small, but the 3D result becomes small 
when Lz is large. 

The reasonable value of may be obtained by comparing the results of 2D mode of 3D 
calculation with calculation in which vertical stratification is taken into account. Tanaka et 
al. (2002) performed such calculations for an inviscid disk. Their method is to decompose 
vertical modes in terms of Hermite polynomials, and they have found that the vertically 
averaged 2D mode gives most of the contribution to the resulting torque. Their 2D mode of 
the potential, V'p^ttwj is given by 

- ™' exp [{x' + y' + 4^) /A] K, [{x' + y' + 4^) /4] , (33) 



2nH 
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where Ko{x) is the modified Bessel function of the zeroth order. Fitting equation (13T1) with 
equation fl33l) . we find = 2.7H gives a reasonable agreement, see Figure [2l Zeroth order 
of modified local approximation incorporated by Tanaka et al. (2002) coincides with the 
shearing-sheet approximation, and we have exactly the same homogeneous equations with 
them for 2D [kz = 0) mode if we assume d/dt = 0. As is clear from Figure [H 2D mode 
in 3D calculation gives smaller amount of torque than the calculation using 2D potential 
given by equation (!29|) . This explains why Tanaka et al. (2002) has obtained slightly smaller 
amount of torque compared to 2D results. We have also checked that setting = 2.7H, 
our results of one-sided torque agree with those determined by using Tanaka et al. (2002) 
methods within an error of 3%. □ 

We use Lz = 2.5H and the mesh number in z-direction is taken to be 128. In to- 
tal, we have performed calculation with {L^, Ly, Lz) = (lOif, 40H, 2.5H) with mesh number 
(512,512,128), and the resulting mesh size is {Ax,Ay,Az) = {0.02H,0.08H,0.02H). How- 
ever, the radial box size is variably extended according to modes in calculation, as discussed 
in Section 12. 1[ Effective box size in x-direction is as much as lO^if. We use eso = 10~^H 
for the softening parameter. Our results vary upto 30% for large values of viscosity when 
smoothing length is varied upto twice the grid resolution. Therefore, our results give at least 
a qualitative view of how disk-planet interaction is altered by the effects of viscosity. The 
variation of one-sided torque as a function of smoothing parameter is further discussed in 
Section 13. 2[ 

3. Results 

In this section, we show our results of density structure and one-sided torque obtained 
by linear analysis. In this section, the normalization of the torque is given by 

TsD = fi'porpH^c'. (34) 

Note the difference in the normalization from that used in Section 12.21 The results of 2D 
mode torque given in Section [3TT] are different by a factor of Lz/ H = 2.5 from those given in 



^ There is another comphcation regarding the normahzation of the torque when comparing the value of 
the torque obtained by our calculation and that given by Tanaka et al. (2002). They use surface density 
to normalize the torque they have obtained. We use normalization given by equation (I32|l . Correspondence 
between the two results is given by reading our pqL^ to Cp of Tanaka et al. (2002), and this is how we have 
obtained the agreement with their calculation. 
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Figure [H The normalization of torque distribution defined by equation fl25|l is taken to be 

T{x,y,z) 



and torque density T(x) defined by equation (!26l) is given by 

T{x) 



(35) 



(36) 



3.1. Calculations Restricted to 2D Mode 

We first show the results of calculations restricted to 2D modes {k^ = 0). Although this 
is only an approximation, physics involved is made clear. Figure [3] shows the torque obtained 
for various viscosity parameter a. For small viscous coefficients {a <^ 10"^), torque is not 
affected by the viscosity, as discussed by Meyer- Vernet and Sicardy (1987). However, when 
viscous coefficient is large, it is shown that one-sided torque increases and peaks around 
a ~ 1. 

Since viscosity, or any form of dissipation, damps density contrast in general, one may 
think that Lindblad torque is a decreasing function of viscosity. [f| However, our result 
shows that one-sided torque peaks at a ~ 1. This result originates from two different effects 
of viscosity. First, viscosity damps spiral density wave, as discussed by Papaloizou and 
Lin (1984) or Takeuchi et al. (1996). The second effect of viscosity may be observed by 
looking at the density structure in the vicinity of the planet. Figure H] shows the contour of 
density profile of the xy-p\ane for a = 10~^ and a = 10^^. It is clear that the oval density 
profile in the vicinity of the planet is slightly tilted for large viscous coefficient, while the 
spiral density wave is damped. This tilt produces asymmetry of density structure in the 
y-direction, thereby exerting one-sided torque on the planet. 

Figure [5] shows the torque density for a = 10^^ and a = 10~^, respectively. It is 
clear that the peak of torque density locates slightly closer to the planet for a = 10~^ than 
a = 10~^. The main contribution to the torque comes from the flow structure in the vicinity 
of the planet, which is evident when we consider the asymmetry of the density perturbation 



The effect of viscosity on the fluid elements trapped in the horseshoe regions enhances the resulting 
torque since viscosity keeps the asymmetry of the potential vorticity. This applies corotation torque, see 
Masset (2001) for detail. In shearing-sheet calculations presented here, however, there is no corotation torque 
since background values are assumed to be constant. 
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between y > and y < regions. In Figure [6], we plot the contour for the sum of torque at 
y > and y < region, 

T{x,y) + T{x,-y), (37) 

as a function of x and y{> 0). This shows the asymmetry of torque distribution in the y- 
directions. In Figure El we compare the results for two different values of viscosity, a = 10~^ 
and a = 10~^. In the case of a = 10~^, the forward-back difference of the torque in the 
vicinity of the planet almost vanishes since the density perturbation around the planet is 
symmetric in the y-direction. In the case of a = 10~^, however, the most of the contri- 
bution to the torque comes from the region very close to the planet since the oval density 
perturbation structure is tilted in y < direction (see also Figure H]) . 



3.2. 3D Calculation and the Effects of Smoothing Length 

In this section, we present the results of 3D calculation. Figure [7| shows the one-sided 
torque as a function of viscosity. The qualitative behavior of torque increasing with viscosity 
is unchanged, but the enhancement of torque becomes large compared to 2D calculation. 
We fit the data and obtain the following empirical relation 

(0.94 + 10a)e-^-^" (38) 



for a < 0.3 and = 2.5H. We have chosen a form of fitting function in such a way that 
torque converges to a non-zero value for a <^ 1, peaks at a ~ 1, and decreases to zero for 
a ^ 1. This fitting formula is in reasonably good agreement with our calculation in this 
range of viscosity coefficient. 

Just as in the calculations restricted to 2D modes, the density structure mainly in the 
vicinity of the planet contributes to the one-sided torque if large values of viscosity is as- 
sumed. Figure [8] shows the torque density profile obtained for different viscosity coefficients. 
The more viscous is the disk, the closer to the planet is the location of the dominant contri- 
bution to the torque. Figure compares the density structure in the yz-plaiae at a; = 0.068if 
for calculations with a = 10"'^ and a = 10"^. The asymmetry in the ^/-direction in strongly 
perturbed region is present in calculations with large viscosity. 

In 3D calculations, gravitational potential in the vicinity of the planet is not as strongly 
softened as in 2D modes. In the vicinity of the planet, gas feels the gravitational potential 
that decreases as — 1/r, where r is the distance from the planet, in full 3D calculations. 
However, the gravitational potential varies only logarithmically when r <^ H (see equation 
fl3T]l ) in calculations restricted to 2D modes. Since the density fluctuation around the planet 
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is given by Sp/po ~ ipp/c^, the deeper gravitational potential gives the higher value of the 
torque exerted on the planet if there is a substantial asymmetry in the density structure. 

Calculations restricted to two-dimensional modes predict one-sided torque very well if 
viscosity parameter is small, since regions close to the planet (r ^ H) is not very important. 
However, if large values of viscosity is used, it is necessary to perform full three-dimensional 
calculations since the density structure close to the planet is important. 

At the end of Section [221 we addressed that varying smoothing length changes one-sided 
torque upto 30%. Below, we discuss the effects of smoothing length in our calculation and 
argue that the values of one-sided torque obtained for large viscosity, a ^ 0.01, seems to be 
the lower limit of the one-sided torque. 

Figure [TO] shows how one-sided torque of three-dimensional calculation varies with 
smoothing length esD. Results with a = 10~^ and a = 10^^ are shown. In this calcula- 
tion, we used the box size with (L^., Ly, L^) = (lOif, lOH, 2.5H) and the grid number with 
[N^, Ny, Nz) = (512,512, 128). Note that resolution in the y-direction is better by a factor 
of four compared to the parameters used in Figure [71 Grid resolution is ~ 0.02H in all 
directions in Figure [TOl 

It is shown that one-sided torque converges well for small values of viscosity parameter. 
If a is as large as 0.1, the calculation converges for sub-grid smoothing lengths, and results 
vary approximately 30% if smoothing length with twice the grid scale (~ O.OAH) is assumed. 
The qualitative behavior of one-sided torque can be understood if we notice that smaller 
smoothing length gives a deeper potential in the vicinity of the planet, and if large values of 
viscosity is assumed, density structure in the vicinity of the planet is important for one-sided 
torque. 

Since contribution to one-sided torque in case of large viscosity mainly comes from 
the density structure close to the planet, it may depend on the grid resolution used in 
the calculation. Figure [11] shows the one-sided torque obtained for different resolutions. 
Calculations with (A^^., Ny, N,) = (512, 512, 128), (256, 256, 64), and (128, 128, 32) are shown 
while keeping (L^., Ly, L^) = (lOiJ, lOH, 2.5H). These values corresponds to grid resolutions 
with Ax = Ay = Az = 0.02H, O.OAH, and 0.08H, respectively. The value of softening 
parameter is kept = 10~^H for all the calculations. It is shown that the value of one- 
sided torque for a = 10~^ is well converged while the value of the torque for a = 10~^ varies 
approximately 30% when grid resolution is varied by a factor of four. We note that in case of 
small viscosity, we have obtained well-converged values of one-sided torque since the effective 
Lindblad resonances, which are located at |a;| ^ {2/3)H, are all well resolved. 

We note that the vertical averaging, large softening length, and coarser grid all introduce 
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more softened potential in the vicinity of the planet and therefore one-sided torque becomes 
smaller in case of large viscosity. Equations fl30|) and fl3T|) show that the vertically averaged 
potential diverges at the location of the planet more mildly. Gravitational force becomes 
weaker in the vicinity of the planet if we use larger values of softening length. The larger 
grid scale cuts off gravitational potential at larger distance. 

From the behaviors when we vary the softening parameter and the grid scale, we con- 
clude that the values of the one-sided torque obtained in our calculation for large viscosity 
is the lower limit. We also argue that the numerical treatment in the vicinity of the planet 
may cause a quantitative difference in the one-sided torque since it may change the density 
structure in the vicinity of the planet. If high values of viscosity is assumed, the density 
structure close to the planet is important and therefore, three-dimensional calculation with 
high resolution and small softening parameter is essential. The resolution and the softening 
parameter should be determined, in principle, by considering the realistic size of the planet. 



4. Analytic Treatment of Density Structure around the Planet 

We have seen that the viscosity exerted on the disk can change the density structure 
in the vicinity of the planet and in a viscous disk, the planet experiences one-sided torque 
from the gas well inside the effective Lindblad resonance. In this section, we show that a 
dissipative force distorts the density structure close to the planet by using a simple analytic 
model. 

We consider a two-dimensional model with a friction force exerted only in the y-direction. 
We consider equations with viscosity terms replaced by friction. We assume k^it) ~ 

since significant excitation of perturbation occurs when radial wavenumber becomes zero. 
We also set kz = and consider 2D perturbation for simplicity. The set of equations we 
solve is ^ 

— -|- iky6vy = 0, (39) 

j^6v^ - 2%6vy = 0, (40) 
5vy + ]^%5v^ = -c^iky + - j5vy, (41) 



and 

d 



dt 

where a = 6p/ po and 7 is a drag coefficient. 

Using equations (139!) and (140|) to ehminate 6vy and integrating once the resulting equa- 
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tion, we obtain 

6v^ = -'^^a + V, (42) 

where is a constant of integration. Physically, V is the perturbation of vortensity, which 
is actually conserved in this model since we neglect terms with kx{t) and we assume friction 
force is exerted only in the ^/-direction. Therefore, V is zero if we assume there is no vortensity 
perturbation initially. Since we are interested in the perturbation that arises from the planet 
potential, we assume V = Q hereafter. 

Using equations fl5^ and fH21) to eliminate dv^ and 5vy from equation fllT]) . we obtain a 
single telegraph equation with a source terms, 

'^^cr da r, ^ 

where 

si = cX + nl. (44) 

The right hand side of equation (l43l) is the source of perturbation that is induced by 
the planet's gravity. The source potential il}{kx{t), ky) depends on time through the time- 
dependence of kx(t). 

Let us now derive the steady state profile of density perturbation in {t, x, |/)-coordinate 
space. In order to solve equation ( l43l) . we Fourier transform perturbation in i-direction 

/oo 
dsf{s,ky)e'-'\ (45) 
-oo 

where / denotes any of perturbation quantities. We note that the Fourier transformation in 
t-direction is equivalent to the inverse Fourier transformation of k^ modes into x-coordinate 
space. The "frequency", s, in equation (HSj) is the frequency of perturbation experienced by 
a single mode specified by ky. 

The real space quantity f{x,y) is obtained by 

/oo roo 
dky dkj{kx,ky)e'^''^'^^''yy\ (46) 
-oo J —oo 

In a steady state, the value of k'^ can be taken arbitrary [see discussion in Section 12. Ij . 
Therefore, we can take k'^ = without loss of generality. Time-dependence of kx is given by 

kx{t) = ^Qpkyt, (47) 

and therefore, we obtain 

3 

dkx = -9,pkydt. (48) 
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Changing the integration variable from to t in equation (H6|) . we obtain 



POO roo 2 

f dky / dt-npkyf{t^ky) exp 

^-oo ^ 

(ifcj^ / dt-^pkyf{t,ky)exp 

-oo J —oo 



2 kyXVtpt -\- kyU 



'kyOC^^pt ~\~ kyXj 



(49) 



The value /(t, ky) appeared in this equation is Fourier transformed in the t-direction accord- 
ing to equation (H5|) . Substituting equation fHOl) into (l45ll . performing integral in t and s, 
we obtain the relationship between /(s, fc^) and /(x, y), 

/■oo 

f{x,y)= 2tt 



dky ^ ^pky 



X 



= -(3/2)npfc^a;, fc,) + e-^^«V(s = {3/2)QpkyX, -ky)] . (50) 



One-sided torque is calculated by equation (but there is no integral in the 2;-direction 
in two-dimensional problem considered in this section). We choose the phase of the potential 
such that ipp{s,ky) is real. It is possible because the potential of the planet is spherically 
symmetric and the planet is fixed at the origin of the coordinate system. Thus, we obtain 



dkyky 



dslm [(t(s, —ky)] 'ippi^s, —ky) 



(51) 



T oc 

Jo Jo 

where constant of proportionality does not depend on ky, 7, and x and "Im" denotes the 
imaginary part. From this equation, torque is exerted on the planet only when a{s, ky) and 
ippi^s, ky) are out of phase. In other words, the amount of torque can be estimated by looking 
at the imaginary part of a{s, ky), provided that the phase of Fourier transform is taken in 
such a way that ipp{s, ky) is real. If density perturbation and the potential are in phase, the 
y-component of the force exerted on the planet is symmetric in the ?/-direction and there is 
no net one-sided torque when integrated over y. 

Now we consider the solution of equation (143|) . By Fourier transform in t-direction, we 
obtain 



a{s, ky) 



IS'J 



klipp{s,k 



y)i 



(52) 



and as equation fl50|) indicates, it is sufficient to consider s = {3/2)QpkyX mode. 

0, a and ipp are out of phase only at the resonance. 



In the limit of small friction, 7 
which is located at 



9q2t2 2 



(53) 



From equation fH4l) . we see that this corresponds to the location of effective Lindblad reso- 
nance (Artymowicz 1993). There is a phase shift of tt only in the vicinity of the resonance 
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and torque is localized at the resonance location. This localization of the torque comes from 
our assumption of zero radial pressure term in equation (HOj) . We note that in this case, the 
resulting torque is independent of the amount of friction if integrated over the resonance 
width (Meyer- Vernet and Sicardy 1987). This can be seen by considering the imaginary part 
of the density fluctuation a limits to, for sufficiently small 7, 

lim Im[a(s, ky)] = lim ^-^—-^I—^ fcj^/^p = 7r5z?(s^ - s^)kliP^, (54) 

where 5d{x) is Dirac's delta function. The last expression is independent of the amount of 
friction, and the torque exerted by the mode ky is determined by the strength of the source 
ip-p at the resonance. Using equations and ( l53l) . we see that there is no resonance close 
to the planet, or region \x\ < Xc, where 

2 c 

In this region, density perturbation and source term are in phase. In terms of real space 
coordinate {x,y), this corresponds to ipp{x,y) oc a{x,y). 

In the case of finite friction, the situation is different. From equation fl52p . we see that 
resonance width is given by |s — so| ~ 7. Therefore, when 7 ~ Qp, the region in which density 
perturbation and source term are out of phase can overlap the location of the planet. Note 
that So is a function of ky and x [see equation (1551) ]. so "the width of the resonance" refers to 
the width in the x-direction, and this width varies with the mode in the ^/-direction (ky) we 
consider. The finite width of the resonance causes the asymmetry in density perturbation 
in y > and y < regions even at the location of the planet, x ~ 0. Since gravitational 
force exerted by the planet on the disk is large in the vicinity of the planet, the significant 
amount of one-sided torque can be exerted on the planet. 

The amplitude of perturbation is suppressed when significant friction is exerted. From 
equation fl52l) . we see that, in the vicinity of the planet 

klijp{s,ky). (56) 



The suppression is significant only when 7 exceeds flp. 

One-sided torque is calculated by equation (but there is no integral in the 2;-direction 
in two-dimensional problem considered in this section). Using equations (15T!) and (!52|) . it is 
possible to obtain 

Toc^ dkyk^yj^ ds-^^^-^^^^^l{s,-ky), (57) 
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where constant of proportionality does not depend on ky, 7, and x. The value of the torque 
is determined by the competition between the suppression of the amount of the density 
perturbation and the amplification of the amount of the torque as the resonance width 
becomes wider. Therefore, it is expected that the torque peaks at 7 ~ Qp. This explains 
the peak of the torque for a ~ 1 in the calculation of viscous disk presented in the previous 
section. 

The qualitative behavior of one-sided torque can be captured more clearly if we consider 
a further simplified case. We consider the case when the forcing potential ipp{s, ky) does not 
depend on s, 

^p(s,y = ^. (58) 

This corresponds to the case where forcing potential is constant in the x-direction. In this 
case, it is possible to perform s integral analytically. 



T oc 



1 



1 



7r + 2 tan 



log 



-1 



2sl 



7 



- 7^ 

7^^M±VZEi3 

72-2sg-7v/7'-4s2 



7' < 
7' > 



(59) 



If there is only one ky mode in the potential, this gives the exact amount of torque. Otherwise, 
this gives the amount of torque exerted by each mode of ky. Limiting values of equation fISIJl) 
are 



T 

T oc 



nk^ 
1 



I7 - 2so| 
T ^ 



7^0 



7 

7 - 



So 



00. 



(60) 



Therefore, the amount of torque is independent of the amount of dissipation in the limit of 
7 — i> 0, it increases as 7 increases and peaks at 7 ~ sq, and then it decreases to zero as 7 
becomes very large. In case of small 7, the integrand of fl57|) is localized at s ~ sq- When 7 
is large, however, it is necessary to consider contribution from all the region of s. In general, 
since the amplitude of forcing term tpp{s, ky) increases as s — > 0, it is expected that the value 
of torque peaks at 7 ~ sq. If the amplitude of forcing term cuts off at ~ ^"^^y^ peak 
of the torque is expected at 7 ~ fip. 
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5. Summary and Discussion 

5.1. Summary of Our Results 

In this paper, we have performed hnear analyses of density structure as a result of grav- 
itational interaction between a planet and a viscous gas disk using shearing-sheet approxi- 
mation. We have obtained the density structure in the vicinity of the planet and discussed 
its effect on type I migration by calculating one-sided torque. We have used time-dependent 
methods to obtain the density structure and a wide range of viscous parameter a is inves- 
tigated. We have calculated one-sided torque, which is the torque exerted on the planet by 
the gas on one side of the planet's orbit. One-sided torque converges to a well-defined value 
if viscosity is small, but increases as we increase the value of viscosity for a ^ 0.01, and then 
decreases as a exceeds the order of unity. The values of one-sided torque we have obtained 
for large viscosity may be the lower limit of the actual value. We note that total torque 
is not actually calculated in this paper since it is not possible within the shearing- sheet 
formalism. However, we also emphasize that the results of shearing-sheet calculations are 
useful in calculating differential torque (Tanaka et al. 2002). We have seen that the torque 
exerted on the planet in a viscous disk is primarily determined by the tilted spherical (or 
spheroidal) density structure in the vicinity of the planet, in contrast to the inviscid case, 
where main contribution of the torque comes from effective Lindblad resonances and the 
horseshoe region. In an inviscid disk, although the value of density perturbation is large in 
the vicinity of the resonance, it is symmetric in y-direction and therefore, torque exerted from 
y > and y < cancels each other. In a viscous disk, there is an asymmetry in y-direction, 
resulting in increasing torque with increasing viscosity until a ~ 1. By using a simplified 
analytic calculations in Section HI we have shown that it is possible to consider these effects 
of viscosity as the increase of the width of Lindblad resonances. When the viscous coefficient 
is larger than a ;^ 1, however, the torque becomes small with increasing viscous coefficient 
since the density perturbation is damped due to large viscosity. The torque can be a factor 
of two larger than inviscid case. The impact of detailed structure around the planet on 
type I migration has not yet been fully analyzed in previous numerical simulations probably 
because relatively large values of viscosity (a ^ 10~^ — 10~^) are required to see this effect. 
Below, we show some discussions and future prospects. 
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5.2. Prospects of Modified Local Analysis and Global Calculation 

5.2.1. Validity of Linear Analysis on One-sided Torque 

In this subsection, we discuss whether or not hnear approximation well describes the 
real density structure around the planet. In the linear approximation, the magnitude of 
perturbed quantities such as Sp/po or 5v/c must be smaller than the order of unity In 
our calculation, all the non-dimensional perturbed quantity is proportional to the value of 
GMp/Hc^, which is about 10^^ — 10^^ for typical values of protoplanetary nebula. Therefore, 
when Sp/p becomes the order of ten in the normalized calculations given in this paper, linear 
approximation becomes less accurate. Since the contribution of the torque in a viscous 
disk mainly comes from such strongly perturbed regions, it may be necessary to perform 
high resolution numerical calculation in order to fully obtain the structure of the density 
fluctuation. This is one of the future prospects of this study. However, we expect that 
the basic physics and the order of magnitude of one-sided torque is similar to the linear 
calculation. 



5.2.2. Differential Lindblad Torque 

The total torque that is exerted on the planet is the differential torque, which is the 
difference of the torque exerted by the inner disk [x < 0) and outer disk [x > 0). In the 
shearing-sheet calculations, the inner and outer disks are symmetric and therefore torque 
exerted on the two regions cancel each other. We have seen that the one-sided torque is 
mainly exerted by the tilted spheroidal density profile around the planet if strong viscosity 
is exerted on the disk. Therefore, we expect that the asymmetry between the inner and 
outer disk comes from the asymmetry of the density structure around the planet, probably 
within, or around, the distance of the order of Bondi radius, tr = GMp/c^. Assuming 
that the differential torque Tdifj mainly comes from the effects of curvature, the order of the 
magnitude is 

Tm ~ Tone-side X , (61) 

rp 

where Tone-side IS onc-sidcd torque obtained in this paper, and 5x is the length scale where 
most contribution to the torque comes from. In the inviscid case, Sx H and therefore 

^^ifT,inviscid ^ Tone—side, inviscid ^ ■ (62) 
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In the viscous case presented in this paper, 5x may be as small as Bondi radius tb- Therefore, 
the lower limit of the differential torque in viscous disk may be estimated as 

-^difT, viscous ~ -^one— side, viscous ^ 5 (6*^) 

and the ratio of differential torque in viscous disk and that in inviscid disk may be given by 

-^difF, viscous ^one— side, viscous '"B f/^A\ 

r T If- ^^^^ 

diff, inviscid one— side, inviscid ^ 

Since 

^ - ^ f^)^ - 8 X 10-3 (J^) (^) ' , (65) 

H Mc \HJ \MJmJ \H/rJ ' ^ ^ 

and Tone-side,viscous/^one-side,inviscid ~ 2 from our calculatiou, differential torque in a viscous 
disk can be much smaller than that in an inviscid disk, although one-sided torque becomes 
larger when the value of viscosity is a ~ 1. 

In order to include the effect of asymmetry, which comes from curvature effect, it is 
necessary to proceed to modified local approximation (Tanaka et al. 2002) or to perform a 
global calculation. In order to find the precise value of the torque exerted on the disk and find 
how various physical processes change the nature of type I migration, the time-dependent 
formulation of modified local approximation is necessary. 



5.2.3. Corotation Torque 

In the shearing- sheet formalism, we cannot calculate corotation torque. For two-dimensional 
modes, this is because all the background quantities (and therefore vortensity) are assumed 
to be constant. For three-dimensional modes, we note that, in Tanaka et al. (2002), there 
is a singularity at corotation modes (see equation (A5) of their paper). However, this singu- 
larity does not account for the torque in itself. Results of shearing-sheet calculations must 
be combined with the higher order solutions of modified local approximation in order to find 
the torque exerted at corotation resonance, see equation (57) of Tanaka et al. (2002). 

If we assume that time derivatives of equations (jH]) and ([7]) are zero in order to obtain 
a stationary solution, we also find a similar singularity at a; = for three-dimensional waves 
{kz 7^ 0), even we have assumed the constant background density in the 2;-direction. However, 
within the shearing-sheet formalism, it does not give a torque at the corotation. Also, this 
singularity does not seem to play a major role in our time-dependent calculations, since there 
is no singularity in our formulation. 
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If we proceed to the modified local approximation, we expect that it is possible to 
calculate the corotation torque in the linear regime. Recently, Paardekooper and Papaloizou 
(2009a) has shown that large viscosity will push the corotation torque into a linear regime. 
Therefore, not only from the point of view of differential torque but also from the view of the 
corotation torque, extension of our time-dependent methods to modified local approximation 
is an interesting future work. 

In summary, in order to predict the precise value of the torque, it is necessary to use a 
modified local approximation or perform global calculations, which can take into account the 
curvature effect. Non-linear analysis may be important since the main contribution comes 
from the region in the vicinity of the planet where density is strongly perturbed. However, 
we believe that the basic physics that exerts torque onto the planet can be captured by this 
linear analysis and therefore, high-resolution study is necessary. In global calculations, we 
also note that it is also possible to investigate the effects of gas accretion onto the central 
star on type I migration, which is always present in a viscous disk but is not captured in our 
local calculations. 



5.3. Flow Structure in the Vicinity of the Planet 

Recently, Paardekooper and Papaloizou (2009b) has shown an interesting results re- 
garding the flow structure close to the planet. In this section, we observe streamline of the 
flow in the vicinity of the planet and compare our results with those previously studied. We 
note that since we have not calculated axisymmetric {ky = 0) modes, the flow structure we 
have obtained is incomplete. However, we can still find some qualitative difference between 
inviscid and viscous calculations. 

Figure [12] shows the streamlines plotted over the density fluctuations for calculations 
corresponding to Figure HI Results with a = 10~^ and a = 10~^ are shown. Calculations are 
restricted to 2D modes of 3D calculation. Mass ratio between the planet and the central star 
q = Mp/Mc and disk aspect ratio h = H/r^ are assumed in such a way that q/h^ = 0.0252, 
just as Paardekooper and Papaloizou (2009b). 

We find that the width of horseshoe orbit is ~ O.IH for a = 10~^, while the width 
becomes slightly narrower in calculations with a = 10~^. We have obtained about a factor 
of two smaller horseshoe width compared to Paardekooper and Papaloizou (2009b) in cal- 
culations with small viscosity (see Figure 3 of their paper.) We think that this is probably 
because we have weaker potential in the vicinity of the planet because of the effective soft- 
ening introduced by vertical averaging, see equation fl^ . In Paardekooper and Papaloizou 
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(2009b), they also obtained narrower horseshoe width in calculations with larger softening 
parameters. However, since we have neglected axisymmetric modes and non-linear effects, 
which may be potentially important in determining the streamline, this issue needs further 
investigation. 

5.4. Axisymmetric Modes 

Although we have neglected the axisymmetric modes in our calculations, a viscous 
overstability in axisymmetric modes has been discussed, especially in the context of Saturn's 
rings. In this subsection, we discuss whether this overstability affects our results. 

Local linear analysis of viscous overstability using a simple hydrodynamic model is 
performed by Schmidt et al. (2001). Stability criterion for isothermal case depends on the 
derivative of viscous coefficient v with respect to surface density perturbation. We note 
that the system of equations we have solved, equations (|2U]) - (125]) . are actually stable against 
viscous overstability, since our assumption of viscosity corresponds to /5 = — 1 of Schmidt et 
al. (2001). In this case, we do not expect viscous overstability to occur and therefore, we 
can safely neglect the axisymmetric modes. 

However, disk may be prone to viscous overstability if different prescription of viscosity 
is taken into account. If the non-linear consequence of viscous overstability produces non- 
axisymmetric structure in the vicinity of the planet, this may exert additional torque onto 
the planet. However, observation of the particle simulation performed by Schmidt et al. 
(2001) (their Figure 1 for example) indicates that non-axisymmetric modes are not strongly 
driven by the viscous overstability and therefore this may not play an important role in plan- 
etary migration. Nonetheless, the effects of viscous overstability seems to be an interesting 
extension of the analysis presented here. 

5.5. Turbulent Disk 

In this subsection, we briefly discuss that our analysis may indicate that the stochastic 
torque is important in a turbulent disk dominated by magneto-rotational instability (see 
e.g.. Nelson and Papaloizou 2004). The effective turbulent viscosity originated in MRI may 
become as much as a ~ 0.1 (e.g., Sano et al. 2004). Therefore, our calculation indicates that 



^ Compare our equations (O and ([7]) and their equation (8). Note there is no self-gravity and bulk 
viscosity in our calculations and we assumed isothermal equation of state. 
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the density structure around the planet is more important than the tidal wave launched at 
the effective Lindblad resonances. Since the density structure around the planet is turbulent 
in the vicinity of the planet, stochastic force may be exerted on the planet embedded in such 
a turbulent disk. 

In order for the a-viscosity prescription [see equation (127|) ] to be a good approximation 
of the turbulent flow, the length scale in consideration must be larger than the eddy size. In 
the present problem, the scale in consideration is smaller than disk thickness. In a turbulent 
flow driven by MRI, eddies of sizes on the order of several tenths of the disk scale lengths 
exist. This can be understood if we notice that the wavelength of the most unstable mode 
is much smaller than the scale height of the disk with weak magnetic field. If weak poloidal 
magnetic field, B^o, is exerted on the unperturbed disk, there always exists a net magnetic 
fiux, 

{B.) = B,o, (66) 

where (Bz) is the horizontally averaged 2;-component of magnetic field. Therefore, there 
always exists eddies with scales of the order of the most unstable mode of magneto-rotational 
instability with net fiux Bzo, 

where A is the eddy scale, vaa = -B^o/V^vrp is Alfven velocity constructed from the average 
poloidal field, and Q is Keplerian angular frequency. In most cases, Alfven velocity is smaller 
than sound speed by more than an order of a magnitude, 

VA,o ^ (68) 

This indicates 

A ^ IQ-^H. (69) 

Thus, we can possibly use the effective turbulent viscosity in order to study the qualitative 
effects of the interaction between a turbulent disk and a planet. 

However, the effective values of a may be obtained by averaging turbulent stress over 
several scale heights, since it includes all sizes of eddy that is important in turbulent fiow 
driven by MRI. Thus, the use of a prescription may not be valid in discussing very small 
scale structure, and full MHD calculation is necessary. More detailed, quantitative analyses 
require high-resolution numerical study of the magnetized turbulent disk. 

We have found that the effects of viscosity becomes apparent if a exceeds the value of 
~ 0.01 — 0.1. Actual values of a in the disk is largely uncertain, but estimated values are 
around this value. Numerical simulations by Sano et al. (2004) indicates that the values of 
a varies from 10~^ — 10^^ depending on the setup. From the observational constraints of 
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dwarf novae, which seem to be better studied than protoplanetary disks, it is indicated that 
the values of a can vary from 0.01 to 0.1 (e.g., Cannizzo et al. 1988). The values of a seems 
to be an open issue for both theoretically and observationally. 
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Fig. 1. — Comparison of the torque obtained in 2D calculation and the 2D mode of 3D 
calculation for various L^. In 2D calculations (dashed line), we use the gravitational po- 
tential given by equation (129|) . while in 2D mode of 3D calculation (solid line), we use the 
gravitational potential given by equation (!3T|) . The value of the normalized torque, T /V2t> 
is plotted, where is defined by equation ( l32l) . 
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Fig. 2. — Comparison of the 2D mode of the potential used by Tanaka et al. (2002) given by 
equation (133|) (sohd hne) and the = mode of the potential given by (13T|) with = 2.7 H 
(dashed line). 
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Fig. 3. — Variation of torque as a function of viscous coefficient a as a result of 3D calculation 
restricted to 2D modes. 
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Fig. 4. — Contour plot of the density fluctuation 5p/ po around the planet (located at the 
origin) for 2D mode of 3D calculation with — 2.5H with a — 10~^ (left) and a — 10~^ 
(right). Axisymmetric mode is not included. 
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Fig. 6. — Forward-back asymmetry of the torque distribution given by equation fl37|) for 
a = 10-^ and a = 10-\ 
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Fig. 7. — Variation of the torque as a function of viscosity parameter a obtained for 3D 
calculations (plus). For comparison, we plot the results restricted to 2D modes by cross. 
Dashed line shows the fitting function given by equation (|38ll . 
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Fig. 8. — Torque density profile obtained by 3D calculation with various a. Solid line shows 
the results of a = 10~^, thick dashed hne a — 10~^, and thick dotted hne a — 10~^. 
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Fig. 9. — The contour plot of 3D mode density structure in yz-plane at x = 0.068H with 
a — 10~^ (sohd hne) and a — 10'^ (dashed hne). The hnes show the contours of 6p/po = 10. 
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Fig. 10. — One-sided torque obtained by 3D calculation for different smoothing length. 
Horizontal axis shows the values of smoothing length e^D and vertical axis shows the one- 
sided torque. Calculations with a — 10~^ (solid line) and 10~^ (dashed line) are shown. 
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Fig. 11. — One-sided torque obtained by 3D calculation for different grid resolutions. Hori- 
zontal axis shows the values of grid size (the same for all three dimensions) and vertical axis 
shows the one-sided torque. Softening parameter = 10~^H is used. Calculations with 
a — 10~^ (solid hne) and 10~^ (dashed hue) are shown. 
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Fig. 12. — Density and streamline for calculations restricted to 2D mode in 3D calculation 
with a = 10~'^ (left) and a = 10^^ (right). This corresponds to the calculations presented 
in Figure m Mass ratio between the central star q = Mp/M^ and disk aspect ratio h = H/vp 
are assumed in such a way that q/h^ = 0.0252, which is the same as Paardekooper and 
Papaloizou (2009b). Gray scale shows the density perturbation Sp/po divided by q/h^ and 
solid lines are streamlines. 



